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Chapter 1 


Introduction 


Safe air is a most vital environmental resource. Its quality is especially critical in an 
enclosed environment such as spacecraft, submarine, airliner cabin. The importance 
of the closely related issue of the indoor (workplace) air quality has been repeatedly 
emphasised in the last twenty years (Berglund et al., 1986; Berglund et al., 1989). It 
is known (Mplhave et a/., 1986) that only a slight degree of air contamination can 
have an acute health and performance impact on humans. 

In addition to background contamination, resulted from material offgassing, mi- 
crobial growth and normal human activities (James and Coleman, 1994), there is 
a risk of an accidental contamination release. In the case of spacecraft and other 
inclosed environments an unexpected event can lead to mission failure or even in- 
capacitation and death since there may be a significant time lag before personnel 
evacuation becomes possible. An early detection of the contamination event can 
provide the critical time necessary to take mitigating actions, or to arrange for safe 
evacuation. 

Recent history of space exploration presents some examples of the possible sources 
of the accidental air contamination (James et al., 1994). One such source is thermod- 
egradation of electronic components in flight hardware (Todd et al., 1993; Todd et al., 
1994). It was reported that during Space Shuttle flight STS-28 0.1 grams of poly- 
tetrafluoroethylene (PTFE, or Teflon by its trademark), used as a wiring isolation 
was pyrolyzed within 1.5 seconds. The crew experienced no adverse health effect. 
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However, it was subsequently estimated that a pyrolysis of 2 grams of PTFE would 
have affected health and performance of the crew. 

Another example of accidental contamination is the fire extinguishant (CO 2 ) re- 
lease due to a containment failure, and microbial growth resulting in the release of 
volatile pollutants into air (Coleman and James, 1994). Based on the experience 
with the Shuttle, it is estimated that minor contamination events will occur with a 
frequency of once in 30 days, and moderate accidents (those causing symptoms in 
the crew) will occur once per year. Details are not available on the accidental re- 
leases into Russian spacecraft. However, it is known that air contamination problems 
have nearly ended a mission (Covault, 1983) and that formaldehyde exposures are a 
concern (Peto, 1981). 

Currently utilized spacecraft air revitalization systems are designed to maintain 
equilibrium concentrations of contaminants within limits prescribed by spacecraft 
maximum allowable concentrations (SMACs) (National Research Council, 1994 and 
1992). The monitoring of the Space Shuttle air is done off-line, and involves sampling 
of the spacecraft air immediately before launch and late in the mission (James et al . , 
1994), with additional samples taken if crew suspects a contamination problem, and 
subsequent analysis after spacecraft landing. The analysis of taken samples shows 
that air quality is consistently within the target region. At the same time, this 
monitoring technique failed to reflect some know contamination problems (such as 
abovementioned PTFE thermodegradation), or provided inconclusive results. It is 
also insensitive to the temporal and spatial variations in cont aminan t concentrations. 
At the same time, it is known that temporal variation of some compounds, such 
as carbon dioxide, can be quite large because the efficiency of the removal system 
degrades quickly when the reactant bed is almost consumed. The spatial variation 
of the contaminant concentration to our knowledge has not been estimated, but one 
can predicted that it can be significant due to nonuniformity of the air flow within 
habitat. This assumption is supported by information (Savina, unpublished data) on 
large spatial variations of the contaminant concentrations throughout the Mir Station 
under normal operating conditions. 
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NASA plans for extended human presence in Space Station and bases on the Moon 
and Mars in the 21st century (White and Lujan, 1989) raises serious questions about 
the adequacy of current air revitalization and monitoring system for long duration 
missions. Added to the risk of accidental contaminations, there is a concern (Logan, 
1989) that Space Station would become an example of runaway tight (or sick) building 
syndrome (Godish, 1995; Mplhave, 1989), since there will not be a luxury of frequent 
returns to Earth for thorough cleaning, characteristic to the Space Shuttle missions. 

It becomes apparent that a new generation air revitalization system should sup- 
port such advanced functions as real-time air quality monitoring; health risk evalua- 
tion of the current situation and prediction of the future risks based on the extrap- 
olation of the current trends; an early detection of a contamination event which in 
turn allows for preventive mitigation of the health hazard; and isolation of the con- 
tamination source, which would facilitate repair and clean-up. Essential components 
of such system are gas sensors, high performance on-board computer and algorithms 
and methods which provide required system functionality. The analysis of existing gas 
sensors technology (Kocache, 1994; Barry and Brackbenbury, 1991) indicate that the 
instrumentation needed for on-line multicomponent gas analysis is already available. 
In fact, a prototype combustion products analyzer has been flown as part of NASA 
program to develop capability to detect combustion products (James and Coleman, 
1994). More advanced instrumentation can be expected to become available in the 
near future, including real-time particle detector suitable for space flights (NASA, 
1989). At the same time, a computing power of available workstations suitable for 
on-board implementations is constantly increasing. 

In this report we present our progress, results and findings in development of 
an integrated air quality modeling, monitoring, fault detection and isolation system 
(Figure 1). The focus of this year efforts was on development of distributed model of 
the air contaminants transport, the study of air quality monitoring techniques based 
on the model of transport process and on-line contaminant concentration measure- 
ments, and sensor placement. In the next chapter we discuss different approaches to 
the modeling of the spacecraft air contamination, and propose a three-dimensional 


3 


distributed parameter air contaminant dispersion model, applicable to both laminar 
and turbulent transport. We further proposed a two-dimensional approximation of a 
full scale transport model based on the spatial averaging of three-dimensional model 
over least important space coordinate. A computer implementation of the transport 
model is considered in the third chapter where we present a detailed development of 
two- and three-dimensional models illustrated by contaminant transport simulation 
results. In the fourth chapter we suggest the use of a well established Kalman filter- 
ing approach as a method for generating on-line contaminant concentration estimates 
based on both real time measurements and the model of contaminant transport pro- 
cess. We also show that high computational requirements of the traditional Kalman 
filter can render difficult its real-time implementation for high-dimensional transport 
model and propose a novel Implicit Kalman filtering algorithm which is shown to lead 
to an order of magnitude faster computer implementation in the case of air quality 
monitoring. Section five discusses sensor placement. In Conclusions we summarize 
our results, and outline the direction for next year’s research. 
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Mitigating Actions 



Figure 1.1: Integrated air quality monitoring system 




Chapter 2 


Dispersion Model of the Air 
Contaminants 

2.1 Lumped and distributed models 

Simulations of a spread, introduction and removal of airborne contaminants creates 
the basis for the solution of the entire complex of problems of monitoring and envi- 
ronmental control of a space habitat. 

Generally, there are two different approaches to simulation. In the first one, each 
separate volume of a habitat, such as a cabin or a room, is represented as a well 
mixed tank with uniform distribution of all parameters. Mathematical models of this 
approach are a system of ordinary differential equations, resulting from the application 
of the macroscopic mass balances (National Academy of Sciences, 1981). The result 
of simulation is time change of volume-averaged concentration of different gaseous 
and airborne particulate contaminants, which depend on the rate of contaminants 
generation and removal. These models are relatively simple to develop and apply but 
are characterized by low resolution and accuracy. 

An example of a computer model of this type is Trace Contaminant Control 
Simulation (TCCS) computer program (Perry, 1993), which was specifically designed 
to facilitate the development of a spacecraft contamination control hardware and 
prediction of its performance. 
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Another example is Computer Aided System Engineering and Analysis (CASE/A) 
modeling package for environmental control and life support system (ECLSS Integra- 
tion Analysis, 1990), developed at McDonnell Douglas Space Systems Company. Al- 
though it also uses a well-stirred tank as an underlying representation of the habitat’s 
volume, it belongs to a class of multi-tank models of indoor air quality (Ozkaynak 
et al., 1982; Ryan et al., 1988), and provides means for simulation of a number of 
connected tanks. This allows to model of different configurations with more than one 
interconnected volumes in habitat or to represent a single volume as a combination 
of well stirred-tanks thus increasing resolution and reliability of simulation results. 

An alternative approach is to base the computer model on the fundamental dis- 
tributed laws of physics. This results in the mathematical description of the process 
in the form of partial differential equations with appropriate initial and boundary 
conditions, and gives an exhaustive information on the distribution of contaminants 
throughout the spacecraft. However, its complexity and high computational require- 
ments are obvious shortcomings. 

In the analysis of the environmental risk of a space mission, the assessment of the 
effects of chronic inhalation of low concentration contaminants is at least as important 
as the analysis of extreme concentrations as a result of on-board accidents. Models 
of the spacecraft modules, based on the assumption of well mixed volumes cannot 
account for fluctuations in concentrations within a habitat. At the same time, effect 
of a higher than average concentration of contaminants and spatial distributing of the 
concentration can be significant, especially for long duration space missions. A dis- 
tributed model is also needed to precisely identify unknown source of contamination. 
All these factors lead to a conclusion that the perspective of long-term and remote 
space missions requires the development of sophisticated models of the space habitat 
which can be run in real time, and serve as a basis for monitoring and decision mak- 
ing systems far more sophisticated than current detection and environmental control 
systems. 
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2.2 Three-dimensional model of the contaminant 
dispersion 


According to Fick’s law of molecular diffusion (Bird et al., 1960), binary mass trans- 
port process with isotropic diffusion can be described by the following convection- 
diffusion equation: 

^ + {VqU) = ( V-D u Vq)+F , 

U = ( u,v,w ), (2.1) 


with an appropriate initial and boundary conditions, where q is the instantaneous 
contaminant concentration distribution within spatial domain ft, Dm is the molecular 
diffusivity of contaminant in the air, U is the velocity distribution of the air flow, F 
is a source function describing generation and removal of the contaminant, and V is 
the vector differential operator. 

If the contaminant fraction in the air is low, equation (2.1) can also be used to 
describe a multicomponent transport. In this case, the spread of each contamination 
component is described by a model in form of equation (2.1). The transport models for 
different contaminants can be coupled through a generation term F, and a molecular 
diffusivity D M i in general depends on a particular type of the contaminant species. 


If the density of the air is constant, equation (2.1) in rectangular coordinates yields 


dq duq dvq dwq 

dt + dx dy dz 

x 


(V • Dm^q) + F, 
{x, y, z) € D. 


( 2 . 2 ) 


An appropriate boundary conditions for equation (2.2) correspond to impermeable 
wall, air duct, open hatch and completely or partially permeable wall (Table 1). For all 
practical purposes, source function F can be adequately described by the combination 
of pointwise contaminant sources and sinks: 

f=E firm* - <o - E tfm* - *?), (2.3) 

1=1 t=l 

where (N so , F*°, *®°) and (N ai , F* 1 , ***) are the number, capacity and location of a 
point source and sink, x = ( x , y, z), and S is the Dirac delta function. 
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Equation (2.2) is directly applicable to the laminar transport of an airborne con- 
taminant. Though results of Son and Barker, 1993 suggest that space station air 
flows will be mainly laminar, turbulent transport can play a significant role at least 
in some parts of the spacecraft habitat. Thus, it is important to modify (2.2) in order 
to accommodate the case of the turbulent transport. 

For the case of turbulent flow, both the flow velocity U = ( u,v,w ) and the 
concentration q must be treated as stochastic quantities. Following the framework 
of semi-empirical theory of turbulent diffusion (Bird et al., 1960) we introduce the 
following notation: 

U = U + U' = (u,v,w) + (u', v',w'), 
q = « + «', (2-4) 

where U, q are time smoothed values and U', q' are instantaneous fluctuations. Since 
we are interested in a meaningful trend of contaminant concentration, rather than its 
stochastic fluctuations, we average the Fickian model (2.2) over a time interval AT 
long enough for the integral of instantaneous fluctuations to vanish. This yields 

^ + V • {W) + V • (qU) = (V • D M Vq) + F, (2.5) 

at 

where an overbar denotes a time averaged value. According to semi-empirical theory 
of turbulent transport, the fluctuation moment q'U' can be approximated by the 
gradient relationship 

tfU 7 = -D T Vq, (2.6) 

where Dt is an empirical coefficient of eddy or turbulent diffusivity. Since 

D m < D t , (2.7) 

molecular diffusion can be ignored, and resulting three-dimensional turbulent trans- 
port model becomes 

^ + V • (qU) = (V • D T Vq) + F. (2.8) 

dt 

It is important to note that the models for laminar (2.2) and turbulent (2.8) 
transport are of the same mathematical form. This allows us to use the same computer 
model for laminar, turbulent and mixed contaminant transport cases. 
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2.3 Two-dimensional approximation 


In an attempt to develop a model which gives sufficient details about contaminant 
concentration profile and at the same time is simple enough to be run on-line, a two- 
dimensional approximation of the transport model (2.2) was previously suggested 
(Todd et al ., 1993; Todd et al., 1994). The transport model can be simplified by 
averaging equation (2.2) over least important spatial coordinate. Assuming the aver- 
aging over z, define the height averaged concentration: 

1 r H 

92 = H Jo qdZ ' ( 2 - 9 ) 


where H is the habitat’s height. Further assume that Dm does not change with z, 
and that (it, v ) is a z-averaged air velocity vector. Integrating (2.2) with respect to 
z, and accounting for boundary condition w(0) = w(H) = 0 obtain 


dq 2 duq 2 dvq 2 _d dq 2 , d dq 2 _ , 
dt dx dy dx M dx dy M dy ’ 


( 2 . 10 ) 


where 


rH 

Jo 


Fdz + Dm 


dq 

dz 


\H 


dq 

dz 


( 2 . 11 ) 


oJ 


Similarly, a turbulent dispersion is approximated by the following two dimensional 


model: 

dq 2 duq 2 dvq 2 _d dq 2 d dq 2 

dt dx dy dx T dx dy T dy ’ 

where q 2 is height averages and time smoothed contaminant concentration. 


( 2 . 12 ) 


This approach provides the flexibility of using different two-dimensional approxi- 
mations (resulted from the averaging along different spatial coordinates), depending 
on the current needs, or to simultaneously run more than one approximation to obtain 
more detailed information about the spread of the contaminant. 
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Chapter 3 


Computer Implementation of the 
Transport Model 

3.1 Two-dimensional model 

The transport equation with an appropriate boundary conditions is used to develop 
a computer model of the transport process. 

First, using finite difference or finite element method find an approximation of 
the spatial partial derivatives in the model equation. For two-dimensional model 
(equation (2.10) or (2.12)) let matrix A is a discrete analog of the spatial operator 

d d d d du dv 
x dx M dx + dy M dy dx dy ’ 

obtained using finite differences. A semi-discrete analog of (2.10) in a vector form 
can now be written as 

ft = Aq+f ' (31) 

where q(t) = [<7i, • • • ,</n] is a finite difference approximation of the concentration 
q 2 (x, y, t ), f(t) is an approximation of the source function f(x, y, t) plus a contribution 
from the boundary conditions. If Cx is approximated using five point stencil, the 
resulted operator A is n-by-n sparse matrix with five or less non-zero elements on 
each row. Dimension n of the system (3.1) depends on the spatial mesh used to 
approximate Cx- 
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Temporal discretization of equation (3.1) concludes the development of the dis- 
crete analog of the transport equation (2.10). In order to ensure numerical stability of 
the computer implementation, for a large range of time steps At it is advisable to use 
implicit approximation of the time derivative in (3.1). For instance, Crank-Nicolson 
approximation scheme results in the following unconditionally stable discrete model: 

Mq m +\ = + fm, (3.2) 

where m is a time index, f m = f(mAt), and 



where I is a unit matrix. Note that if coefficients of (2.10) are time dependent, 
matrices A\ and A? will change with m. 

The sparsity of the system (3.2) can be exploited to yield an efficient implemen- 
tation of the computer model. The use of the indexed storage of matrices A\ and A 2 
leads to a significant reduction in the required computer memory, while utilization of 
the structure specific methods of matrix manipulation results in an efficient solution 
engine for linear system (3.2). 

The details on development and implementation of two-dimensional transport 
model based on 5-point stencil approximation of spatial derivatives, Crank-Nicolson 
approximation in time domain, and spatial solution methods, which take advantage 
of a sparse structure of the discrete transport model are given in the next subsection 

The finite elements method leads to the similar results. A spatial discretization 
in this case yields the following semi-discrete equation 

da 

M^ = A fe q + f fe , (3.3) 

where M and Aj e are usually called mass and stiffness matrices. They are sparse, 
and depend on weight and element functions of the chosen finite element method, 
and on the ordering of the components of q. A source term f j e is an approximation 
of f(x, y, t) plus a contribution from the boundary conditions. 
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As before, a discrete analog is obtained as a result of the temporal approximation 
of the above semi-discrete equation. An implicit discretization yields a discrete model 
in the form of equation (3.2), where A x and A 2 depend on the specific scheme. One 
possible expression is 


A 1 = 


A 2 — 


(M 
\At 
( M 
\A t 



obtained using central difference approximation in time domain. 


3.1.1 Computer implementation of two-dimensional trans- 
port model 

The model equation (2.10) of the two-dimensional transport with appropriate bound- 
ary conditions is used to develop a computer model of the transport process. First, 
using finite difference approximations of the partial derivatives with respect to spatial 
coordinates we obtain a semi-discrete analog of the model equation. The diffusion 
operator is approximated as 


dx 2 


^nj 


A* 


<?n-H,p 2(?n,p 4~ Qn- l,p 

Ax 2 

9n,p+l ^9n,p ~b <?n,p— 1 


(3.4) 

(3.5) 


dy 2 *'" ,p A y 2 

where q n<p ~ q 2 (t, nAx,pAy), Ax and Ay are discretization steps along coordinates x 

and y , V n<p corresponds to V = {D M (nAx,pAy),D T {nAx,pAy)}, and the subscript 


is used to specify a point on the spatial mesh {(n,p) | n = l,N,p = 1,P}. 

Since central difference approximation of the convective operator is known to cause 
non-physical oscillation of the numerical solution, the following upwind differencing 
scheme is used: 

_ ,,w 

(3.6) 

OX - L\x 

where 


duq 2 

. 9n,P - C 

dx 

~ Ax ’ 

,E _ , 

f 9n,p ^ ^n,p ^ 

n,p 

\ Qn- l,p if ^n,p ^ 0? 
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w _ 1 

Qn,p j 

[ 9n— l,p ^ ^n,p ^ 0> 
l Qnj, if «n,p < 0. 


Similarly, 

dvqi 

dy 

<p~<v 
V n ’ P Ay ’ 

(3.7) 

where 

«n,p = | 

[ 9n t p if v n,p ^ 0, 

[ 9n,p— 1 if ^n,p ^ 0? 



II 

Cr 

f Qnj>—1 if Vn,p > 0, 
[ 9n,p if ^n,p ^ 0- 



After spatial discretization the model equation take the following semi-discrete 
form: 


dt 


a E -q w 
Vn,p Vn,p 


q s — 


+ “~ a v 


_ -r» ( Qn+l,p %Qnj> + Qn-l,p , Qnjrt-1 ~ %Qnj> + Qn,p-l\ , em 

~ n ’ p \ Ax 2 + Ay 2 J + W 


(3.8) 


The application of center difference approximation of the time derivative yields 
the following discrete analog of the two-dimensional transport model: 


q m + 1 -q m 1 

Vp q JhB. = _± , u. 


At 


n,p 


r „e i m+1 

r„wi m+1 

\ n E l m 

r Wl m 

[? n ,p 

L^.pj 

\- 7/- 

l^n.pj 

l v n.pj 


Ax 


Ax 


+ v, 


® n,p 


r ? 5 i m+i _ u i 

Ay 


m+1 


+ v n, } 


k,r- ter 

Ay 


' m+l _ 9 „m+l . „m+l -m+1 _ o n m+l , n m+ 1 
9n+l,p Z( ?n,p + Vw— l,p + 9nj>+l ^g n ,p + «?n,p-l 


Ax 2 


+ 


C+i,p - + C-1.P | <g+i - 2^p + Cp-1 


Ax 2 


Ay 2 


Ay 2 

+ c. 


(3.9) 


where At is the time discretization step, and the superscript m = 0, 1, 2, ... is used to 
identify an instance t—(m+ l)At for which the solution of equation (3.9) is sought. 

After collecting like terms in (3.9) obtain the following equation for a single spatial 
mesh point (n,p): 


A<c*!„ + &C‘ + CCi + D,“, + = 
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(3.10) 


~Mn-l ,p + G Qn,P ~ C <ln+l,p ~ D C, P -l ~ 1 + fn,p, 

where the coefficients are equal 


A 

= A x 

+ a 2 , 

B 

= B x 

+ B 2 , 

C 

— 2?n+l,p7 

D 

= A 

+ e> 2 , 

E 

— £*n,p+lr 

G 

= A 

+ g 2 , 


Ay 


Ax 


where 


and 


1 Ay 

A ' = 

AxAy Vn+t^pAy 2> nj>+ iAy 
1 At 2Ax 2Ax 


( Ay Ax \ 


Di = 

G\ = 


V n>p Ax 


2Ay ’ 

AxAy _ Vn+x^Ay _ Z> n , P+ iAy 
At 2Ax 2Ax 


2Ax 

_ v + 

l 2 Arc 2 Ay 


A 2 = 

d 2 = 
b 2 = 
b 3 = 


l h 

g 2 


_t*n l£ Ay < n 

2 ? u n,p ^ 

Un .x p.Qi u >o 

2 i u n,p ^ u > 

Vft .P Ag < q 

2 5 u riyp ^ u ) 

Vn iF Ag 7; >0 

x 2 5 Un # ^ U ’ 

B% + B\ 

u » .p A y n <o 

2 j u n,p ^ u ) 

_ u ",p A y u >o 

2 j u, n,p ^ 

yu? A .£ <jt <0 

2 i y n,p ^ 

-^xP.^g 7; >0 

x 2 ’ V *»P ^ U ’ 

= G3 + G4, 


(3.11) 


(3.12) 
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(3.13) 


Ga — 


■{ 


1 !a ^,u w >0, 

v "<p Ax v >0 
2 j u n t p 

The discrete transport model (3.10)-(3.13) is unconditionally stable and is free 
of non-physical oscillations. It has a second order accuracy in the time domain. 
Though, application of upwind differencing for convective operator results in only 
first order accuracy in space, the absence of numerical oscillation for all values of 
model parameters and discretization steps makes model (3.10)-(3.13) very attractive 
from the practical point of view. 


3.1.2 Discrete model in the matrix form 

The system of N x P equations of the form (3.10)-(3.13) describes a contaminant 
transport process within a spatial domain 0 < x < (N — l)Ax, 0 < y < (P — l)At/. 
However, its exact appearance depends on how these equations are arranged into a 
single matrix equation. We have implemented a node numbering algorithm known as 
D4, or altering diagonals method (Aziz and Settari, 1979), which allow a particularly 
efficient direct solution of the model equation. According to this method, an equation 
describing each node (n,p) enters as i-th row of the overall matrix equation (3.2), 
where determination of i is summarized in Table 1. 

The discrete transport model can now be written in the matrix form as 

Aiq m+ i = b , (3.14) 

where vector q m+1 = {q™* 1 | 1, N, 1, P} approximates the contaminant concentration 
distribution on the time step m + 1, A x is a sparse NP x NP matrix with at most 
five non-zero elements on each row, and 

b=A 2 q m + f m , (3.15) 

where f m accounts for boundary conditions, and sources and sinks of contamination. 
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Even (n + p ) 

Odd (n + p) 

For (n + p) < N 
. = (2±£-l) 2 +p 


For N < n + p < P 

i - ULU 2 
1 4 

(N— l) a — JV(n+i>)— 1 
2 

For N + l<n + p<P + i 

■ _ NP N 2 

1 ~ 2 4 

+ n±£ N + l-n 

For n+p> P 

■ _ (N- 2) 2 . NP-N 2 -2N 
‘ — 4 1 2 

. (2N+P-n-p+4)(n+p-P) 

For n + p > P + 1 
i = NP — P + p 

(N+P— (n+p)+l)(iV +P—(n+p ) — 1) 

1 4 

-(n + 1) 

4 


Table 3.1: D4 Node Numbering Algorithm 


3.1.3 Solution of the matrix equation (3.14) 


The alternating diagonals node numbering method produces the discrete model (3.14) 
which can be partitioned as 


-^11 

A12 


" 9 m+l " 


hi' 

A21 

A22 


Qm+ 1 . 


.^ 2 . 


(3.16) 


where An, A 22 are diagonal matrices, and A 12 , A 21 are sparse matrices. 

Special structure of equation (3.16) suggests a simple and efficient method of direct 
solution of matrix equation(3.14). First, by direct exclusion transform the lower half 
of the matrix equation (3.16) to obtain 


'An 

A 12 


9 m + 1 


hi' 

0 

A22 _ 


Qm+l. 




(3.17) 


A sparse matrix equation 

A 22Qm+i = b 2 


(3.18) 


of reduced order &£■ s now independently solved for q^+i using Gaussian elimination. 
Given the solution of equation (3.18), 9m + i is found as 


Qm+i ~ ■4-n^i A n Ai2<7 m+ i, 


(3.19) 
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where An is a diagonal matrix. Overall solution of (3.14) is now formed as 


Qm+l 


" 9m+l " 
,9m+ 1 . 


(3.20) 


3.1.4 Simulation example. 

For the cabin depicted in Figure 2 consider a two-dimensional height averaged ap- 
proximation of the contaminant transport process. Assuming the laminar air flow, 
the parabolic input and output air velocities, and additional parameters presented in 
Table 2 we can calculate 2 - averaged air velocity distribution in the cabin as shown 
in Figure 3. Here arrows follow the streamlines of the air flow, and there length is 
proportional to the magnitude of the air velocity at a particular spatial point. The 
maximum air velocity is equal 1 m/s and occurs in the center of the air duct. Note 
that due to a large velocity gradients we have assumed V = 23 cm 2 /s, which is the 
eddy diffusivity of carbon dioxide in air. 



12.2 m 


Figure 3.1: Simulated cabin geometry 

Let us further assume that initially contaminant concentration in the cabin is zero, 
and that at the time t = 0a contaminant is introduced into cabin with the inlet air 
stream at a known rate (Table 2). 

Using this contamination scenario, we can calculate the evolution of the contami- 
nant concentration distributions in the cabin. Results of the simulation 6 sec, 40 sec, 
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12.2 m 


Figure 3.2: Height averaged air velocity distribution 

2 min, 4 min, 6 min and 10 min 40 sec after the beginning of the emission are 
shown in Figure 4. At 10 min 40 sec contaminant distribution has already reached 
its steady state. 

3.2 Three-dimensional model 

Described two-dimensional design can be used to discretize the partial differential 
equation (2.1) or (2.8) to obtain a straight forward implementation of the three- 
dimensional transport model. However, the dimension n of the resulted system of 
the form (3.2) can be very large. At the same time, the structure of the sparse 
matrices A\ and A 2 is more complicated in this case. Specifically, the width of the 
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band in Ai and A 2 can be substantial. Consequently, strait forward implementation 
three-dimensional transport model can be quite computer intensive. 

A significant reduction in computational complexity can be achieved following the 
paradigm of the alternating-direction implicit (ADI) approach, which embodies the 
idea of operator or time splitting. Such splitting allows the reduction of the original 
problem into a series of simplified problems which must be solved consecutively on 
each time step. For example, the Douglas method leads to the following numerical 
implementation: 

A x — Q = {^A-x + + 2 A z — Qm "F 

(~ A » - J») 9 " = - it q ‘’ 

(321) 

where q* and q** are intermediate variables, Ai is a finite difference approximation 
of the spatial operator C>, : i = {x, y, z), where, for instance, 

, . d d du 

1 dx M dx dx 

Thus, instead of solving one large and complex system (3.2) we need to solve three 
tridiagonal equation (3.21), which can be accomplished simply and efficiently. 
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Chapter 4 


Air Quality Monitoring Using 
Kalman Filtering 


A modeling of the introduction, dispersion, and removal of airborne contaminants in 
an enclosed environment is very important. It summarizes our knowledge about the 
transport process, and is especially valuable during the design stage, since it allows us 
to predict what happens under different scenarios of chronic or accidental contamina- 
tion, analyze the influence of different factors on the spread of contaminants, simulate 
performance and limitations of air revitalization system, and use pre-recorded flight 
data to recreate and analyze real events occurring on-board of the spacecraft. At 
the same time, the model is merely a reflection of our current a priori knowledge 
about physics of the air contaminants transport process, geometry of the volume of 
interest, emission sources and removal devices. As a result, the quality of simulation 
depends solely on the assumptions employed in the development of the model, and 
completeness and accuracy of the input data. 

In reality we almost never have all the required information and the adequateness 
of the model is often in question. Furthermore, a real system is subject to stochastic 
perturbation and information about some system parameters (such as eddy diffusiv- 
ity) is inherently uncertain and incomplete. All this decreases the value of the results 
based solely on a priori knowledge and input data. 

The most important information about air quality in a habitat is the concentration 
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distribution of the airborne contaminant, pre-selected for monitoring. The task of 
generating a real (or near real) time contaminant concentration estimates can only 
be accomplished using on-line concentration measurements. However, one should be 
aware of the limitations of the measurement data along. Practical considerations allow 
to deploy only a limited number of sensors. Each sensor produces spatially localized 
(i.e. pointwise) information, which is corrupted by noise, and requires filtering. The 
measurements from different sensors may not be in complete agreement requiring data 
fusion and reconciliation. And finally, localized (pointwise) measurement data must 
be extrapolated over the entire habitated volume to obtain complete contaminant 
concentration distribution. 

All these serves as a motivation for developing the air quality monitoring system 
based on both on-line measurements, and the verified model of the transport process. 
The integration of the measurement data and the model can be achieved within the 
framework of the well established Kalman filtering method. According to Kalman 
filtering paradigm the uncertainties of the model and measurements are represented 
by the additive stochastic white Gaussian perturbation. First, lets consider the mod- 
ification of the transport model. For the transport model in rectangular coordinates 
(2.2), the modified stochastic model of the process takes the following form: 

^ + duq + dvq + dwq = {v , D M V 9 ) + p + ^ t ) w ( aSj *), (4.1) 

at ox ay oz 

?(®»0) = q Q {x), 

E[g°(x)] = 0, E[g°(x 1 )g°(* 2 )] = Pq{xu * 2 ), 

X,X\,X 2 e Cl, 

with an appropriate boundary conditions, where c(x, t) is a deterministic function of 
noise intensity, w(x,t) is Gaussian white (in time) process with zero mean, and 


E[tw(®i,t)tu(* 2 ,i r )] = Q{xi,x 2 ,t)S(t - r), 

where E[-] and £(•) are the expectation operator and the Dirac delta function, Q and 
p 0 are nonnegative function, symmetric in the sense that Q(xi,x 2 ,t) = Q(x 2 ,x\,t), 
and p 0 (xi,x 2 ) =p Q (x 2 ,x 1 ). 
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The process model (4.1) must be further supplemented with the model of the 
measurement system, reflecting our knowledge on how the actually available sensors’ 
output relates to the state of the process q. In a general form, measurement system 
can be represented as 


z(x,t) = h(q,t) + v(x,t), x € f), (4.2) 

where z(x, t ) is the output of the measurement system, h(q, t) is a spatial measure- 
ment operator (usually linear), and v(x, t) is Gaussian white process independent of 
w(x , t) with zero mean and nonnegative symmetric covariance function R: 

E[v(xi,t)u(x 2 ,T)] = R{xi,X2,t)6(t - r), xi,x 2 € Q. 

Currently available sensor usually can provide only pointwise readings. A general 
model of the measurement system (4.2) in this case take form 

z{Xj,t) = h(q(Xj, t), t ) + v(xj, t ), j = 1, . . . , l, 

where Xj is a sensor’s location, and l is a total number of sensors. 

The level of the model and the measurements uncertainties is determined by the 
covariance functions Q and R. The choice of the measurement noise covariance matrix 
R is made based on the careful study of the sensing instrumentation in the controlled 
conditions. The covariance of the model noise is chosen at the model validation 
stage and depends on how well the model reflects the real process. After initial 
determination of Q and R they are often viewed as a design parameters, and can be 
further adjusted to obtain desired estimates. 

A theory of the Kalman filtering for distributed parameter systems (Tzafestas, 
1982; Ray and Lainiotis, 1978) can be readily applied to the system (4.1)-(4.2). 
However, resulted distributed Kalman filter equation, and an associated nonlinear 
distributed Ricatti equation must eventually be transformed into a discrete form to 
enable a numerical solution of the filtering problem. 

Alternatively, the Kalman filtering can be applied to a discrete analog of the 
distributed parameter system (4.1)-(4.2) to begin with. Though not as theoretically 
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sound, this approach presents some practical advantages since it is less mathematically 
involved, and is more general in the sense that it leads to an implementation of the 
Kalman filter which is less dependent on the changes in the process model. 

The discretization of the stochastic model (4.1) concludes in a discrete analog in 
the form of either a single implicit equation (3.2), or a system of implicit equations, 
such as the system (3.21). 

First, consider the case of implicit equation (3.2), appropriately modified and 
supplemented with the discrete analog of the measurement equation (4.2): 

Ai q m+1 = A 2 q m + f m + C{rn)w m , (4.3) 

Qo = 9 °> 

Zm + 1 = H(m + l)<Z m+ i + v m+ i (4.4) 

where C{m) and H(m+l)q m+1 are discrete approximations of c(x,m At) and h(q, (m + l)At), 
and z m+ 1 is a measurement vector corresponding to z(x,t ); w m and u m +i are uncor- 
related zero mean white Gaussian sequences such that 

E [wjWn] = Q(m)6(j - m), 

E[u j+1 v£ +1 ] = R{m + l)<5(j - m), 

where Q{m) and R(m + 1) approximate functions Q(x i, x 2 , mAt) and R(x i, x 2 , ( m + 1) At), 
and 5(j — m) = 1 if j = m and zero otherwise. 

The traditional Kalman filtering require that the model be given in the explicit 
form. If matrix Ai in non-singular for all m, the implicit model (4.3) can be rewritten 
in an equivalent explicit form, and the traditional Kalman filter can be applied to 
generate the optimal estimates of q m+1 . However, there is a strong motivation to 
avoid matrix inversion step in a filtering algorithm. If for some m matrix A\ is 
ill-conditioned, its inverse is calculated with significant error, unless some special 
measures are built into the filtering algorithm. This would usually involve on-line 
calculation of the condition number of A x and application of iterative improvement 
in order to calculate the explicit model equation with acceptable accuracy, or the 
use of a singular value decomposition as a first step in calculating matrix inversion. 
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Furthermore, the system matrix A^ l A 2 of the equivalent explicit representation must 
be treated as general (full) matrix, since inversion destroys matrix sparsity. Therefore, 
application of the traditional Kalman filtering can leads to significant computational 
errors should the matrix A\ be ill-conditioned, and results in an inefficient algorithm 
due to the need to manipulate full (non-sparse) matrices. 

We have previously proposed an implicit Kalman filtering method (Skliar and 
Ramirez, 1995a, b), which does not require matrix inversion, and when applied to a 
sparse implicit system is an order of magnitude faster than the traditional Kalman 
filter, making it a superb method for on-line estimation of the contaminant concen- 
tration. 

According to the implicit Kalman filter method, the optimal concentration esti- 
mate q m+ i| m+ i is given by the following implicit Kalman filter equation: 

Vm+\\m “k Ty (771 -f- 1) H\{rn ■+• l)y m +i|m] (4-5) 

for m — 0,1, ..., where 

Vm+\\m f m (4-6) 

with q 0 | 0 = q°, and the modified the measurement matrix H\{m + 1) is determined 
by the linear equation 

H 1 A 1 = H. (4.7) 

The implicit Kalman filter gain L y (m + 1) is equal 
L y (m + 1) = P v m+l \ m+1 Hl (m + 1) JT^m + 1) 

= n + .|„X("* + + 1) + «(m + 1)]:}4.8) 

where the predicted estimation error covariance matrix P v m+l \ m of an auxiliary vari- 
able y(m+ 1) = Aiq m+l is found as a result of the time propagation of the estimation 
error covariance P q m \ m of the concentration q m according to equation 

n+n™ = A 2 Pl |ra Af + C(m)Q(m)C T (m). (4.9) 

The estimation error covariance matrix P^ +1 , +1 satisfies the covariance measure- 
ment update equation 

n + i lm+ , = n + „ m - n + ,im»r <»> + 1> 
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x[H,(m + l)PS, +1 | m Hf (m + 1) + R~'(m + l)]~‘H,(m + l)f* +ll i4.10) 
or equivalently 

n+iH+i = V-L,{m+ 1 )H,(ra + 1)] PJ, +1|m . (4.11) 

The error covariance matrix -P^, +1 | m+1 of the concentration estimate g m+1 | m+1 is 
related to the covariance matrix Pm+i\ m +i by the following linear matrix equation: 

n +1 , m+ . = (4.12) 

The implicit Kalman filter generates the minimal variance estimations of </ m+1 , 
and is theoretically equivalent to the traditional Kalman filtering, provided the inverse 
of matrix Ai exist for all k. However, it provides a basis for a new implementation 
algorithm for the implicit system (4.3) which does not require matrix inversion. This 
makes it a superior approach when matrix A\ is ill-conditioned or sparse. 

We now formulate an algorithm to determine the optimal estimate of the state 
q m+1 of the system (4.3), (4.4) using the implicit Kalman filter. Given z m+ 1, q m | m 
and L y {m 4- 1), 

1. Compute y m+1 | m by propagating q m according to equation (4.6). 

2. Solve the linear matrix equation (4.7) for the modified measurement matrix 
H i(m + l). 

3. Solve the linear equation 

-^l^m+llm+l = 2/m+l|m L y (k + 1) ^ m +i _ H\ (m + l)l/ m +i| m ] (4-13) 
for the optimal estimate 9 m+ i| m+ i- 

If matrices A\ and H are time invariant, matrix H i needs to be calculated only once. 
The Kalman gain X y (m+ 1) can be calculated in the following ways. Given P q m (m , 

1. Compute P v m+l \ m by propagating P q m \ m according to equation (4.9). 

2. Compute the implicit Kalman filter gain using equation (4.8). 
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3. Calculate Pm+i\m+i according to equation (4.11). 


Note that in order to initiate the gain calculation algorithm on the next time step, 
we need to find -P/[ +1 |*+i from the linear equation (4.12), using direct or iterative 
methods of solution. 

Turning attention to a discrete analog of the stochastic transport model (4.1)-(4.2) 
in the form of the system of implicit equations we immediately observe that after 
appropriate modifications system (3.21) can be written a single implicit equation 


where 


and 


AiQ m+ i — A-iQ m + 

2/m 

0 

+ 

2C(m) 

0 


0 


0 


Wm, 



("A, " A) 


= {A?} = 


A 2 = {A*'} = 


_ 2 _ 

At 

0 


0 0 

(-A. - s) 0 


2 _ 

At 


0 0 (A x + 2A„ + 2A* - £) 

0 0 ~Ay 

0 0 -A z 


(4.14) 


z m + i = [0 0 H(m + 1)] Q m+1 + v m+ i. (4.15) 


The system (4.14)-(4.15) is in the same form as (4.3)-(4.4), and the implicit 
Kalman filter is directly applicable. Resulted algorithm can be further simplified 
if a special structure of (4.14) is taken into account. After obvious transformations, 
the estimation of the contaminant concentration is determined from the sequential 
solution of the following tridiagonal equations: 

(“A* ” 9* = (A* + 2 A„ + 2A 2 - — ^ q TO(m + 2 f m + L x [z - i?iy m+1 j TO ], 
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- A >-i i)<”— *■* 


mm 


—q** + L 2 [z-H iy m+1|m ], 


( 2 \ 2 

— A 2 — 9m+l|m+l = ~ J ‘^zQm\m ~ ^9 + Lj[z — H^iy m +i| m ]i 

where the predicted estimation of the auxiliary variable y = A X Q is equal 


( 4 . 16 ) 



4 13 

“2 


2/m 

2/m+l|m 

A? 

9m|m “l” 

0 


Af 


0 


( 4 . 17 ) 


The modified measurement matrix H\ = {H\j} is found from the following linear 
equation: 

[Ha H n H 13 ]A: l = [0 0 H], 

solution of which reduces to the solution of the following three linear tridiagonal 
equations : 

' H u Af = H, 

H u Af = -H n Af , 

HnAl 1 = -H l2 A?. 

The implicit Kalman gain L m+ 1 = [L[ L \ Lj] T is equal 

+«( m + d]' 1 ’ 


( 4 . 18 ) 


(4.19) 


where 
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m+l|m> 


( 4 . 20 ) 


( 4 . 21 ) 


and 


py — 4 , P^ 4 t 


-{[■ 

estimation error of Q. 


where P£ +1 , m+1 


pQ 

* m+l|m+l 


J* J | , i f j = {1,2,3} is the covariance matrix of the 
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Taking the advantage of the block-matrix structure of A\, the solution of the last 
equation is reduced to a sequential solution of the following six tridiagonal equations: 


4'^ 

A! 1 [P« 


a; 1 

Af 


m+l|m+lj 

pQ 

** m+l|m+l 


] U 4 1T = 

pv 

_ m+l|m+l 

11 

) 

]V T = 

pV 

_ m+l|m+l 

12 

: 13 Af T = 

pv 

_ m+l|m+l 

13 

:22 a 22T = 1 

pv 

1 m+l|f7»+l 

1 22 


-A 2l l P® 

A \ 1/ m+l|m+l 


l u 4 2lT 

m-fl|m+lj "1 > 

,Q ] 12 A32T 

m-f l|m+lj 

>e 

m+l|m+l 


1 u a? t 


AT [P m +l|m+lJ 

Af[p; 


Q 1 23 4 33^ _ \pv l 23 - A 21 \p Q ] 1Z A 3 

A l - ["* m+l|m+lj A 1 [■* m+l|m+lj A \ 


)"at t 
- [aT [p 


Q 

m+l|m+lj 
12 -32 T 


i 12 a 22T 


-A 22 \ P Q 1 A 32T - A 21 

A 1 [■* m+l|m+lj A \ A 1 [•* r? 


Q 

m+l|m+l 


Q 

m+l|m+lj 

l 33 *33? Ipy - A 32 \P Q I 22 *32? 

J m+l|m+ lj [** m+l|m+lj 

l 23 4?3 r _ \aP\i>Q . I 23 


1 13 A 33T 


j32fnQ 1 433^ 4 3 2 T pQ I A 3 

"l [■* m+l|m+lj A 1 A 1 [■* m+l|m+lj A 1 


33 ? 


(4.22) 


where [P£ +1 , m+1 ] 33 = ^ +1 | m+r 

We are now in the position to formulate an algorithm for estimation of the contam- 
inant concentration q m+i based on the measurement data (4.15) and the transport 
model (4.14). Given z m+ u q m )m and I'm - f 1 j 


1. Compute y m+1 | m by propagating q m \ m according to equation (4.17). 


2. Successively solve three tridiagonal matrix equations (4.18) for the modified 
measurement matrix H\. 


3. Successively solve three tridiagonal equations (4.16) for the optimal estimation 

9m+l|m+l • 

The calculation of the gain L m + 1 of the implicit Kalman filter follows the following 
algorithm: 

1. Calculate P y m + X \ m according to equation (4.20). 

2. Calculate the implicit Kalman filter gain using equation (4.19). 

3. Calculate P^ +1 , m+1 according to equation (4.21). 
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4. To initiate the gain calculation algorithm on the next time step, sequentially 
solve six tridiagonal matrix equations (4.22) for the error covariance matrix 

pQ 

* m+l|m+l * 

Despite the complex appearance, presented implicit filtering algorithm is very 
effective for the large dimension dim q m = n, and is convenient for computer imple- 
mentation. If the number of sensors l < n, then the computer implementation of 
each the implicit filtering algorithm requires on the order of n 2 floating point opera- 
tions, while the traditional explicit filter will require 0(n 3 ). Additionally, all linear 
systems that need to be solves in the course of implicit filtering have the same simple 
tridiagonal form, and a single solution engine can be used to implement the proposed 
algorithm. 

Finally, it is well known, that for real-time applications of the Kalman filtering 
it is advisable to use a square root implementation algorithm, which has a superior 
numerical stability. The details on the square root implementation of the implicit 
Kalman filter can be found in Skliar and Ramirez (1996). 
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Chapter 5 


Sensor Placement for 
Contamination Detection 

A discrete time and discrete location algorithm for state estimation and the placement of 
monitoring devices (sensors) in a distributed parameter contaminant simulation is described, 
implemented, tested in a simple contaminant scenario, and compared to Extended Kalman 
filtering results. The state estimation portion of the algorithm is a suboptimal variation of 
Extended Kalman Filtering. The sensor placement portion of the algorithm is based on the 
minimization of the trace of the prediction error covariance matrix. It utilizes variables 
from the suboptimal state estimation portion of the algorithm with a yes/no type of 
switching algorithm to optimize the placement and number of the sensors in a set of pre- 
specified locations. 


5.1 Introduction to sensor placement 

Determination of the optimal number, type, and placement of sensors in a habitat is 
important to both ground-based applications and future space applications. "Tight" 
terrestrial buildings do not allow the free exchange of building air with the outside air, 
resulting in a buildup of stale air and potentially harmful toxicants in the building. 
Monitoring devices in the building can alert the building supervisors to high levels of 
toxicants and indicate when the building must be ventilated, but these monitoring devices 
must be placed in locations which can yield appropriate and useful information. It is 
presently not feasible from either a cost or a logistics stand-point to place sensors 
everywhere in a building that designers or operations personnel may deem them desirable. 
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Space habitats can be considered "tight" buildings, but with the added constraint that 
they cannot be ventilated with relatively clean outside air, the time-honored mitigation 
measure for ground-based buildings. It is therefore important to detect potentially harmful 
buildups of toxicants earlier than in terrestrial buildings, and alternative mitigation steps 
must then be taken. The increased detection requirements in space habitats increase the 
complexity of the monitoring system; but once again, sensors cannot be placed everywhere 
in the habitat. In addition to the obvious cost of the sensors, the cost of the information 
processing equipment, the limits on the physical placement of sensors in the habitat, and the 
cost to transport the mass associated with the sensors (from the Earth to the space habitat) 
motivates the space habitat designer to minimize the number of sensors used in the 
monitoring system. 


5.2 Korbicz sensor location algorithm basics 

A prospective algorithm for sensor location has been developed. It is based upon the 
algorithm, developed by Prof. Jozef Korbicz (1986, 1988, 1991) and, uses Kalman Filter 
system observation techniques to place a limited number of sensors at selected locations in a 
two-dimensional area. This analysis technique can be expanded to three-dimensions but 
with more intensive computations. (The contaminant concentration in a compartment 
volume has been traditionally averaged over the height of the compartment so that a two- 
dimensional analysis can be performed on a three-dimensional space.) The algorithm 
minimizes the contaminant concentration prediction error at selected points in the analysis 
area, based on the sensor readings of a specified sensor configuration, and thereby increases 
the estimation accuracy of the contaminant concentration. The algorithm is able to 
determine a near-optimal placement of a specified number of sensors in specified locations, 
given the dispersion of a contaminant as described by a distributed parameter model, and 
given the desired concentration monitoring accuracy. The suboptimal algorithm of Korbicz 
is reportedly more computationally efficient than previous algorithms due to its 
approximation of the prediction error covariance matrix. 
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A distributed parameter model is a mathematical model of a system that is described by 
a set of partial-differential equations. This is contrasted to a lumped parameter model, 
which is a mathematical model of a system that is generally described by a set of ordinary 
differential equations. The distributed parameter model attempts to describe the system in 
terms of the contaminant dispersion and removal processes in both time and spatial location, 
hence the dependent variables are functions of more than one independent variable, and 
partial differential equations are required. The lumped parameter model attempts to 
describe the system in coarser terms of the bulk flows and movements of contaminants in 
the system as a function only of time, discounting the distances between and within the 
"well-mixed" lumped nodes; hence, the dependent variables are only functions of time, and 
only ordinary differential equations are required. 

While the Kalman Filtering algorithm can be used on both lumped parameter systems 
and distributed parameter systems, the filtering algorithm described below is developed 
around and will use a distributed parameter system. But in order to simulate the spread of 
contaminants in the model of the space habitat at discrete times, the system model has been 
basically reduced to small lumped nodes. The concentrations at these nodes are functions 
of both time and the distances to the adjacent nodes. 

The equations have been coded for the suboptimal filtering algorithm (J. Korbicz, 
1986), and for the sensor placement algorithm described in subsequent papers (J. Korbicz et 
al., 1988; J. Korbicz, 1991). The suboptimal filtering algorithm differs from a standard 
Kalman filter in that it approximates the prediction error covariance matrix through a 
discrete, time-stepping algorithm. The suboptimal filter equations generate state estimates, 
prediction and filtering errors, and the prediction error covariance estimator from given 
measurements. 

The sensor location algorithm determines the optimal placement of N’ spatially-fixed 
sensors to be distributed over N potential locations (N*N', maximum of one sensor per 
location) in a 2-dimensional spatial configuration to monitor the distribution of one chemical 
species concentration. The methodology of the measurement location optimization is to 
minimize a cost functional, based on the trace of the prediction error covariance matrix. 
The sensor location algorithm uses data generated by the suboptimal filtering algorithm to 
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generate adjoint variables, which are the dynamic constraints for the minimization of a cost 
functional. 

The sensor location algoiithm uses a yes/no decision algorithm for the optimal 
placement of N' sensors in N locations. A switching function, based on the Hamiltonian of 
the optimal sensor problem, determines whether a sensor is or is not placed at a specified 
location. With the switching function having yes or no as the decision options, this 
measurement problem can be considered as a "bang-bang" control problem. 

Inputs to the Korbicz Sensor Location Algorithm include the possible sensor locations 
within the habitat, the habitat system’s distributed parameter equations (contaminant 
generation and spread), the level of noise that a sensor generates, the level of "process 
noise" from the habitat, and the required monitoring accuracy. Outputs from the Korbicz 
Sensor Location Algorithm are the minimum number of sensors required to meet a specified 
filtering accuracy (one aspect of a minimum cost sensor configuration), and the specific 
placement of those sensors to achieve that accuracy. 

The system process noise from the habitat is the amount of contaminant concentration 
variability acknowledged in the dynamics (generation and spread) of the contaminant in the 
habitat due to airflow variability, disturbances, etc. 

Example 

A simple two-dimensional example will be used to demonstrate how an optimal sensor 
location algorithm works. Assume that the habitat atmosphere layout is rectangular as 
pictured in Figure la. The atmospheric system is gridded to allow observation of 
contaminant concentration over spatial locations within its whole area at n =30 locations 
(grid intersections). A contaminant generation and dispersion model is applied to the 
habitat. The possible sensor locations are the five circles ( N = 5). These are all of the 
permissible locations of sensors, denoted by the location vector x* = 
((1,3),(2,1),(3,5),(4,3),(5,2)). The locations may have been chosen for a variety of 
engineering and other reasons. 
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Figure la 

i = (l 0 1 0 0\H = S 



Figure lb 

A = (0 1 0 1 0\H = 6 



Figure lc 

i = (0 1 0 1 l}# = 3 


Figure 1 . Movement of Sensors Using Sensor Location Algorithm 

The algorithm continues by guessing at the number and placement of sensors in a 
sensor configuration. The habitat model has a certain amount of confidence, reflected in the 
system process noise, and each sensor has a defined noise level (noise generated by the 
sensor). Our initial guess is to place two sensors. A' 7 = 2, at the locations marked with X's 
in Figure la, corresponding to the vector A = (l 0 1 00). The placed sensors are used to 
monitor the habitat in a simulation, and an estimate of the concentration levels at selected 
points in the habitat (states) are calculated. The adjoint variables are calculated and the 
accuracy with which the sensors monitor the contaminant concentration can be judged by 
comparing the prediction error to a pre-determined error threshold. 

Other monitoring parameters can be used as auxiliary measures, such as comparing the 
estimated contaminant concentration levels against the known levels of the contaminant 
concentration at the points, calculating the percentage of maximum peak concentration 
predicted versus the known percentage of maximum peak concentration, or calculating the 
Hamiltonian of the system. For our initial placement of sensors, the scalar Hamiltonian, H , 
is calculated to be 8. Maximizing or minimizing the Hamiltonian of a system corresponds to 
minimizing or maximizing a cost functional (such as the prediction error), according to the 
discrete maximum principle of Pontryagin. This will be described in more detail later; for 
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now let us assume that we want to minimize the prediction error, and therefore minimize 
the Hamiltonian. 

Using a cost function based on the state prediction errors, the sensors are then shifted 
(via yes/no choices) to other locations, A = (0 1 0 1 0), to provide a possible increase in 
accuracy until shifting of the sensors no longer helps, as shown in Figure lb. The resulting 
error of this "best" configuration is then tested. The resulting value of the Hamiltonian in 
our example has been reduced to H = 6 . If the desired accuracy (comparing predicted or 
estimated concentration states against the model's actual concentrations) is attained using 
some computed configuration with that number of sensors, then the number of sensors is 
decreased and the algorithm rerun until the desired accuracy is unattainable. If the desired 
accuracy is unattainable with using only two sensors (our desired accuracy threshold 
corresponds to H < 4), then three sensors can be tried, A = (0 1 0 1 1), as shown in Figure 
lc. This testing and retesting is continued until one of the previous successful 
configurations is revisited. The resulting set of sensors is a feasible suboptimal 
configuration of the number and locations of the sensors as determined by the Korbicz 
algorithm, and the algorithm stops. Our example found the Hamiltonian to be H = 3 for our 
third configuration, which reflected a minimized error in the state prediction estimate. 


5.3 Simplified space habitat model 

A simple, distributed-parameter, contaminant transport system model based on a model in 
one of Korbicz's early papersJ. Korbicz, 1986 was used as a baseline model. The simplified 
model consists of a tubular space, with eleven equal-distant "grid" points at which the 
contaminant concentration levels are desired to be known or estimated. This is illustrated in 
Figure 2. The points are Ax = 0. 1 distance units apart. This is a very simple model, but can 
be considered to be similar to the interior of a cylindrical habitat, inside of which a 
contaminant is dispersing. In this model, the primary transport mechanism considered is 
diffusion, with a small factor of removal of the contaminant. This small contaminant 
removal might be considered as deposition of the contaminant on the walls of the habitat. 


37 


X 1 x 2 x 3 x 4 x 5 ><6 x 7 x 8 x 9 x 10 X 11 


Figure 2. Simplified Discrete Habitat Model 

The transport of the contaminant within this non-linear system is described by the 
continuous model. 


frkO -D ^ykt) b yfet) , . , , 

at ai J i+IjM ^ ' ^ 

(change) = (diffusion) - (deposition) + (control) + (noise) (5 j) 

where y(x,t) is the contaminant concentration vector at locations x at time t. The discrete- 
time model of the continuous system of Eqn. (1) can be approximated by, 

y (* , * + 1 ) = N ^ (y, A: ) + k ) x «(*, k ) + B(s, k)xw($,k) , (5 2) 


where. 


N , (y ,*, A) = y (jf, A) + A/ x 


| 0x £^*) 

1 fa 


-Ax 


j<j,A) y 
1+lKsf.^J’ 


A($,k) = AtxI, 
BQ[,k)=I. 


(5.3) 

(5.4) 

(5.5) 


The system was initialized at concentration values ofy(x,A = 0) = 0.4. The diffusion 
coefficient was set to D = 1 and the removal coefficient was set to b = 0.05, essentially 
removing the removal term from affecting the contaminant concentration (but still 
accounting for it mathematically). The discrete, white, gaussian system process noise, 
wQcJc), was specified to have a mean of zero, 


E[w(x, k)] = 0. 


(5.6) 
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and a (discrete) variance of <J^ = 0.0016, resulting in a covariance matrix of. 


Q($,xl,k) 


0.0016 

0 


0 

0.0016 


(5.7) 


The values of Q reflect the fact that the standard deviation of the model uncertainty of 
each state is o w = 0.04 , or 4% of the steady-state values of the states, which are at 

concentrations ofy(x,k = <*>) = ~1 at steady-state, wfe,*) is also independent over time. 

A step input contaminant flow to the system was input at the xi l location at time t = 
0.01 sec., and the contaminant concentration was held at unity for the duration of the 
simulation, for, 


J&11»*^°-°1)“1- (5.8) 

Other than the step input, no "control" inputs were used (u(x,k) = 0). The equations used 
to calculate the linearized response are (A is a tri-diagonal system matrix), 

y(x> t) = A x y, (5.9) 


where. 



0 


-LL 

A * 2 


0 


(d22-h.\ -XL 
W 2 / A * 2 


(5.10) 
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The removal term. 


bx 


jfek) 

i+|y(i>k)r 


from Eqn. (3) was linearized into b x 


b(x,k) 

2 


by letting 


ly(x,£)| = 1, the idealized steady-state value of ally. The last row of A was set to zeros only 
for the simulation of the step input in the contamination scenario so as to not diminish the 
step input level; the non-zero system model values were used for the input to the Kalman 
Filter. 

Because this simple system is so well connected, it is observable using only one sensor 
at any location. This was tested by examining the rank of the observability matrix (G.J. 
Smith, 199S). This system brings out an important aspect of testing for observability. It 
was necessary to normalize the system matrix of Eqn. (10) to accurately test for 
observability; otherwise, the rank of the observability matrix was calculated to be less than 
full because of the cutoff tolerance of the rank testing subroutine. Also, even though this 
system is observable with using only one sensor, that does not mean that all of the system 
states can be accurately estimated, as we shall see later in this paper. 


5.4 Extended Kalman filter algorithm 

The Kalman Filtering technique is a good method to obtain state information about states 
which are not directly measured, and is used as a baseline algorithm for the example. The 
Extended Kalman Filter is used for nonlinear systems, as compared to a basic Kalman Filter 
which is used for linear systems. Using the methods outlined in (W.F. Ramirez, 1994), an 
Extended Kalman Filter was written, using a least squares estimator, to estimate the 
contaminant concentration state at each of the eleven locations. Candidate sensor locations 
are at x\, *6, XI 1- The sensors were specified to have a noise variance of 0.0025 (which 

yields a standard deviation noise level in the sensor measurements of 0.05, in the units of the 
concentration being measured). 

Using the Extended Kalman Filter algorithm, steady-state values of the prediction error 
covariance matrix and the Kalman gain matrix were calculated, using groupings of two 
sensors or three sensors, distributed one sensor per candidate sensor location maximum. 
Using those steady-state estimation matrices in the filtering equations, estimates of the 
states at each location were then calculated. The estimates of the states were compared to 
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the "actual states", and the absolute values of those differences were summed and averaged 
over the duration of the simulation. This average estimate error magnitude can be 
considered one measure of the accuracy of the Kalman filter. Note that the trace of the 
prediction error covariance matrix, calculated for each configuration, is not a good 
representation of effectiveness of the estimation when one sensor configuration is compared 
against another configuration. The base equations for the Extended Kalman Filter are 
presented below, adapted from 1994). The system and measurement models are given by, 

y(x,k+ 1)= F(y,*)+ B(x,k)xw(x,k), ( 5 . 1 1 ) 

z(*) = N h (f,*)+v(*), (5.12) 

where >>(£,£+ 1 ) is the updated system state vector, z(A) is the measurement vector, wfe,*) is 
the white gaussian system process noise vector, B(k)is a multiplicative factor matrix, and 
v(k) is the white gaussian sensor noise vector. The non-linear system is modeled by the 
vector function F (y,k). The non-linear measurements are modeled by the vector function 
Nh(v,A). The system model variables are assumed to be defined at the locations x, and the 
measurement model variables are assumed to be defined at the locations x*. The Extended 
Kalman Filter equations are. 


y(x,k) = F(p,£,* - 1)+ B(k - l)x w(x,k- 1), ( 5 . 13 ) 


P(k) = M(k)~ M(k 


dy(k) 


M 


cW*) 


\ ( 




dN H (*) 

dy(k) 


\T 

) 


>1 


+m 


'2M& 
, dyi k ) 


N 

M(k), 

J 


(5.14) 


K(k) = P(k) x 


<3N h (*) 


M 


xR~\k), 


(5.15) 
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M(k + 1) = 


, $K*) 


xf(it)x 




M 


+ *(*)x{X*)x* (*)> 


(5.16) 


i(i,*)=Ki^)+^)x(z(*)-N H a*)). 


(5.17) 


The vector y( x, k) is the system state estimates before measurement (based upon the 
previous estimate), and y(x, k) is the system state estimate vector after measurement. K(k) 
is the Kalman Filter gain matrix. P{k) is the state error covariance matrix after 
measurement, while M(k) is the state error covariance matrix before measurement. 


5.5 Extended Kalman filter results 

The accuracy of the Extended Kalman Filter was first tested using sensors at each of the 
three candidate locations, which yielded a time-averaged state error of 0.01343 (in the units 
of the concentration). The state error is the average difference between the actual 
contaminant concentration values, obtained from Eqn. (11), and the estimated state 
concentrations, obtained from Eqn. (17), including the concentrations at all eleven state 
locations in the average. For this sensor configuration, the trace of the prediction error 
covariance matrix, tr(Af), was 0.03256 (in the units of the concentration squared). Table 1 
shows the results using the Extended Kalman Filter. 


Table 1. Results of Extended Kalman Filter 



Sensor 

Average 


Placement 

State 

A = 

\) A k) *<*)] 

Error 

[i 1 1] i 

0.01343 

_ [1 1 0] 1 

004205 

[1 0 1] 

0.01247 

[Q i U 1 
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Next, only two sensors were used to supply the measurement information required by 
the Kalman Filter to estimate the states. The worst estimates were obtained when sensors 
were placed at locations xi & X6, which yielded a time-average state error of 0.04205. 
Sensors were then placed at locations x\ & X] i, and the resulting time-average state error 
was 0.01247, a significant improvement. The best estimates were obtained when sensors 
were placed at locations x6 & xn, which yielded a time-average state error of 0.01199. 
The plotted actual model concentrations and the estimated concentrations for this sensor 
configuration are shown in Figure 3. This sensor configuration was even better than the 
simulation run using three sensors. This is probably because the weightings being used 
place too much emphasis on the sensor measurements in our simulation. Heavier reliance 
on sensor measurements are good during disturbances, but not good during steady-state 
conditions. 

From these results, it is seen that it is important to have a sensor at the location of the 
contaminant source (highest concentration). With one of the two sensors placed at xu, 
placing the other sensor at location x6 yielded only slightly better results than placing the 
other sensor at location xi . 

5.6 Korbicz suboptimal filtering algorithm 

The algorithm for the suboptimal placement of the sensors is in two parts, the suboptimal 
filtering algorithm equations and the suboptimal sensor placement algorithm equations. 

5.6. 1 Objective Function 

The purpose of the suboptimal filtering equations is to calculate the estimated value of the 
states, denoted by y , such that the estimate is very close to the actual value of the states, y. 
This is primarily accomplished by calculating a filter gain matrix, denoted by K, which 
adjusts the computation of the estimated states. The filter gain matrix is based on the 
minimization of the discrete linear Kalman filter performance objective function, J, 

^ x at 1 x{> {& - j*-' xCz- ] (5)8) 
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actual non-linear and estimated (- -) states (xl. x6. xll) 



Figure 3. Extended Kalman Filter Estimate of States, Sensors at *6 & *1 l. 
Response to a Step Input at Location x] i 

where H is the measurement matrix, R is the covariance matrix of the sensor noise, Z is the 
actual measurements, and M is the prediction error covariance matrix (before 
measurement). The first term of this objective function minimizes the deviation between the 
new state and the state estimate before the new measurements are taken, while the second 
term minimizes deviations between the actual measurements and their computed values 
(W.F. Ramirez, 1994). 

From the minimization of Eqn. (18), we obtain a value for M, which leads directly to 
the computation of the filter gain K (that equation is stated later). The product of the filter 
gain matrix K and the innovation error (difference between the actual measurement and the 
predicted measurements) is then added to the system model predicted states to yield new 
values for the estimated states. In short, the filter gain matrix adjusts the system model's 
predicted states based on the measurement errors to yield better state estimates, which are 
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then fed into the system model again. This algorithm, extracted from (J. Korbicz, 1986; J. 
Korbicz et al., 1988; J. Korbicz, 1991) is described in more detail below. 

5.6.2 Algorithm Equations for a General System 

A system can be described by n states within a connected open spatial domain £2 , where the 
states are at spatial locations, x, within the spatial coordinate vector, x, in the spatial domain 
of (as before). For simplicity, we shall assume that the number of states equals the 
number of spatial locations, and therefore consider only one contaminant. 

The states of the system, y (an w-dimensional state vector), are discretized at iteration 
step times k. The time iteration runs from 0 through K. The states are calculated using a 
system model function, N x (state update differential operator vector), and incorporating 
system process noise. The system process noise weighting matrix function is denoted by B, 
and w is the white gaussian system process noise of the individual states (the noise inherent 
in the dynamics of the system). The state of the system at time *+l, y(x,k+l) is equal to 
the update differential operator given information from time k, N x0\£,*), plus system noise, 

Bw 


y(&,k+l)=N x (y,x,k)+ B(js,k)xw(x,k), k = 0,...,K (5.19) 

The covariance of the system process noise is given by Q, and is assumed to be uncorrelated 
across time. 

Out of the n states, N states (discrete measurement points in the coordinate spaces) are 
measured at times k, and the measurements are denoted by z, the observation vector at the 
selected points. The selected physical locations are denoted by x*, which is a subset of x. 
The measurements z incorporate sensor noise denoted by v.. For linear systems, the 
measurement matrix H transforms the states y into measurements z. For non-linear systems, 
this transformation is denoted by Nh(M), a vector function at the measurement points. 

=N hC y,k)+v(k)~Hxy(x,k)+v(k). ^ ^ 
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The state vector at the measurement points is denoted by/*. 


/(*.*) = [>{*'. *)x*2’* • 


(5.21) 


The white gaussian measurement noise v (the noise generated by the sensors) at the 
instrumented locations is derived from R, the specified covariance of the measurement 
noise, and is uncorrelated between sensors. 

The state estimates, y, are calculated from the state propagation function, Nx, and 
adjusted by the product of the filter gain, K, and the measurement prediction error, v (also 
called the innovation error). 


y(x,k+ 1|* + 1)= N x (f,s,*|*)+ K(x,k+ l)x v(k+ 1) . 


(5.22) 


The innovation error, v , is a vector which calculates the difference between the actual state 
measurements, z, and the predicted measurements, N„(y,k+l|k), 


v(* + l)= z(A + l)-N H (y,* + l|*). 


(5.23) 


The initial state estimate at the points in the system is set to y 0 , 

.P(s,A|*L=o =To(*) 


(5.24) 


The system filter gain (weighting) matrix, K, is derived from the product of the prediction 
error covariance matrix, M, the calculated trends of the state measurements as predicted by 
the model, (N H ). (a Jacobian matrix), and the inverse of the T matrix. 


K(x, k + l) = M(y, S ,k + l|*)x (N h ); x r" 1 (k + 1| *), 


where. 
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(5.26) 



dN H i 

dy 1 


<?N H O’ .•**.*) «Nh(m«^)... 

^ n h0’-** ,*) 

#£)>*) #(* 2’*) 

#(*„>*) 


: 




<?n h (^*„.*) 

dyi?\,k) 

#(*«>*) . 


The T matrix includes the covariance of the measurement noise, R, and the prediction error 
covariance matrix, M, adjusted by the predicted state measurement gradients, (N H )- 


r(* + l|*)— R(k + 1) + (N h x M{k + l|i)x (Njj ^ 27 ) 

If the measurement matrix H in Eqn. (20), has only one non-zero numerical value in each 
row, and zeros elsewhere, then each measurement, z, corresponds to an individual state. 
Therefore, the jacobian matrix (N H )» from Eqns. (25) & (27), and expanded in Eqn. (26), 

will have the same values as H. 

The state prediction error, y(x,k + l|k) = y(x,k+l)-y(x,k + l|k), is calculated using 
the state estimate error, y(x, k| k) = y(x, k) - y(x, k| k), updating it by the calculated trends 
of each state, (NJ- (a Jacobian matrix), and including a noise component. 


y($,k+ 1|*) = (N x )^ xy(x,k\k)+ B(x,k)xH>(x,k), 


(5.28) 


where. 



dN„ i 

dy T W 


3N x (h) 3N x fa) 3 N x (m) 
i Nxte) : 

dN x (») dN x (y w ) 


(5.29) 
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The state filtering error, y(x,k+l|k+l), given by Eqn. (30), is comprised of the state 
prediction error, y(x,k + l|k), as adjusted by a filtering term, which is the product of the 
filtering gain, K, and the innovation error, v. 


y(jf,Ar+ 1|£ + 1) = y{&,k+ 1|*)- K(x,k + l)x v{k + 1) . 


(5.30) 


The initial state estimation error, y(x, k|k), is set to be zero. 

.P(£,*Mt=o = 0 . (5.31) 

The one-step-ahead prediction error covariance matrix, M, is approximated, which reduces 
the required calculations. It utilizes the previous calculation of the prediction error 
covariance and adjusts it by a literal interpretation of the error covariance using the product 
of the prediction errors. 


E £v(* +\)-y(k+ 1|*)) (K*+ 1)-X*+ l|*)) T J= e [£(* + l|*)y(* + 1 *) T ]= M ^ 32) 

Since this filtering algorithm converts the system being analyzed from a distributed to a 
lumped-parameter-type model, and thereby reduces the spatial cross-correlation of the 
prediction errors, a "sliding-mean" (to be described later) is utilized to incorporate 
neighboring prediction errors in the calculation of each prediction error. 

M{x,xl,k + tk)= M(x,xl,k\k - l) 

+ Jn + l|*) x 5y T (i', * + lj*) x dCidCl' - A|* - l)j , 

(5.33) 

where Sx and Sx' are the correlation intervals around x and x ’ , respectively, and y(k) = 

The initial value of the prediction error covariance matrix is set to a value somewhat 
dependent on the model being analyzed, though it should be set to a value larger than the 
expected steady-state value. 
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k]/c-\j k=x = Mq (sX) - 


(5.34) 


These equations have been assembled in a MATLAB based code. 

The time averaging factor, y(k) = was chosen specific to the iteration initialization 

of the algorithm code. Since k = 0 on the first pass through the calculation of M from Eqn. 
(33), Y( k )lk-o is initially equal to and the first calculated value of M is equal to the 

average of Mo and the integral term. 

The purpose of the jacobian matrix used in Eqn. (28), and expanded in Eqn. (29), is to 
take the amount of change predicted for the states (contaminant concentrations) over the 
next iteration and apply that proportionally to other variables related to those states, such as 
the state estimation error y(x, k| k), to obtain predicted values of those variables. 

The portion of the state estimation algorithm which sets this method apart from other 
state estimation algorithms is the estimation of the prediction error covariance matrix, as 
given by Eqn. (33). The approximation of M incorporates a "sliding mean" into the 
calculations. A sliding-mean across the space coordinates is used to calculate the individual 
elements within the prediction error covariance matrix so as to partially reintroduce the 
spatial correlation component of the covariance matrix back into the filtering calculations, 
which were removed from the calculations when the system was "lumped" into discrete 
points. The sliding mean used in this case is a spatially moving average of neighboring state 
prediction errors, encompassing spatially adjacent neighboring points from the physical 
habitat being analyzed. The estimated sample prediction error covariances within the range 
of the sliding mean are added together then averaged. 

5.7 Korbicz suboptimal sensor location algorithm 

The following equations for the Suboptimal Sensor Location algorithm take the original 
suboptimal filtering equations and incorporates a switching variable, X, which denotes 
whether or not a possible sensor location in the Euclidian space actually has a sensor, as 
determined from the proposed sensor configuration being investigated at that time (J. 
Korbicz, 1991). X is specified at each possible sensor location, with a value of 1 if a sensor 
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is present, or a value of 0 if no sensor is present in that location. The discrete adjoint 
(costate) variables are then calculated from the associated Hamilton function. The filtering 
and adjoint variables are used to calculate a switching function, which determines the new 
values of X (new sensor locations). 

The revised filtering equations are approximately (only presenting those equations 
which have changed and neglecting some unnecessary subscripts), 

z(*,*) = [*(*,, j = 1,2, 

(5.35) 

The new measurement vector is the same as before, except that each possible sensor 
location (from the N locations) lacking a sensor (at time k) has a value of A, equal to zero, so 
no signal comes from that location. Mathematically, the constraint that only N’ sensors can 
be placed in the ^locations can be stated as, 


X* 


A M =0orl 


(5.36) 


For our space habitat example, we have added the constraint that a maximum of one 
sensor can be placed at any possible location. 

The prior adjustment to the state estimate formed by the product of the filter gain 
matrix and the innovation error in Eqns. (22) and (30) have been transformed into more 
complex equations due to the addition of the sensor switching variable for each possible 
location, A . 


y (f,k + 1|* + 1) = N x (y,*,*|*)+ G(x,k + 1|*), 
y (* , * + ik + 1) = fa * + tf )- G M + 1|* ), 
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(5.39) 


G(s,* + l|*)= + xA a+1 xr '(pjjffik+Tfc) 

i=\j=\ 

X [z£ j,k + 1)- X jM] X N h (y„s y , k + 1| *)] 

The variables are, 

Ajjc = 1, if a sensor is located at the point xj at time k; 0, otherwise. 

N' = number of sensors. 

N = number of possible sensor locations (N 3 N'). 

Nh = vector function at the measurement points, non-linear function of states. 

(N h )* = gradients of state updates (Jacobian matrix). 

The scalar performance functional, J (in continuous form), to be minimized for the 
sensor location portion of the algorithm is the integral of the cost function, 3, which has 
been selected to be the trace of the prediction error covariance matrix, 

J = l t0 tr[Ar(i,a.',0}* ( 540 ) 

This minimization seeks the reduction of the prediction error over the period of time being 
analyzed. The trace of the prediction error covariance matrix is the sum of the squares of 
each state's prediction error. This captures the largest values in the matrix, and yields a 
positive value for each. Minimizing the trace reduces the error between the predicted states 
and the actual states, and places the sensors at the locations which are associated with the 
least overall system state estimation error. 

The scalar Hamilton function of the system, H (in continuous form), is the sum of the 
cost function, 3, and the product of the adjoint variables, <|>, t|t, 0 , times the right hand side 
of the suboptimal filtering equations (the discrete forms are given by Eqns. (22), (30), and 
(33)). The Hamiltonian in continuous form is, 
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(5.41) 


+ t{e\ x x,-)x^f^\ 


The corresponding adjoint variable equations as calculated by Koibicz are. 


N' N' 


= -(N x J x<t>(x,k + \\k+\)+^'£(N u ty,x J ,k + lk'fy. xA >.*+ 1 

«= 1 7=1 ; 

x r _1 (jf y ,$,, * + l|*)x A a+1 x (N H (y ,*„* + l|*))- 
x M t (x, x f ,k + H^)x \y(x,k + l|fc+ 1)- <p(x,k + 1|£+ 1)], 


V(x,k\k) = — — 


M3 

m 


-3^5: x 9 (i,i,* + |*+i)xj n HgMyteL, 


W,1,J 3Ar() 


= -/ + y(k) xO(x,j i',k+ 1 ) 


XX £ £7 (* + 1 )>* + 1 )- ^ 7 ,*+i xN h (P,jl # 7 (k + 1 ),*+ 1 |*)J 


1 = 1 7=1 


(5.42) 


(5.43) 


xr l (sj (k + 1 ),£.',. (k + 1 ),* + |*)x A, 

,*+1 x (^h(f^,(* + + 

x[VC*,*+ 1)- <t>(x,k + l)]x/ , 


(5.44) 


with initial conditions. 


0(i> k )k=K ~ ® » 
V(.X,k) k=K =0, 


(5.45) 

(5.46) 
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6(x,iL,k) k=K = Q, 


(5.47) 


where, 

<j>, \| r, 0 = adjoint variables, 

tr[*]= trace of a matrix (sum of the main diagonal), and 
K = final iteration count. 

The above adjoint variable equations are integrated backwards from the ending time k — K 
to k = 1, calculating the switching function for each point. The switching function as 
calculated by Korbicz (J. Korbicz, 1991), is derived from the portion of the Hamiltonian of 
Eqn. (67) that is linear with respect to X. The switching function as calculated by Korbicz is 
as follows, 

e jk =e[[ (x, k + \\k + 1)- y» T (x, k + \\k + 1)] 

x M* ^Xj(k+l),k + l|*)x (n h (k,<,(*+ i),*+i|*)) r xr 1 £,(*),*+ 1|*) 
x[z^(* +1 U+l)-NH^£:y(^+lU+lk)]l> = 1 - 2 ’-- ,X,k=l,2,... .,K 

(5.48) 

The switching fiinction values, Ejk, are averaged over time (due to the E[«] term) for 
each point. The summed values are then ordered from maximum to minimum, in the order 
of the best estimated locations for the sensors. The N locations which have the largest 
values of Ej are then considered to be the best locations given the restrictions on the number 
of allowable sensors (for that run). These are the sensor locations to be set (Xj = 1) in the 
next iteration run; sensors will be placed in those locations. 

The optimal sensor locations are updated after each run until no further improvements 
are apparent, as judged from the trace of the prediction error covariance matrix. It is then 
decided whether the optimized sensor configuration will yield sufficiently accurate results 
for the contaminant scenario. Korbicz states that the trace of the prediction error 
covariance matrix can be used for this. Another measure of the accuracy of the sensor 
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configuration is the average error between the estimated states and the actual states. This 
average state error measure was used by the author in this study, as it gives a better 
estimate of the errors that will be produced by the filtering algorithm. 

If the optimized sensor configuration does not yield sufficiently accurate results, then 
the number of sensors used in the configuration is increased after each run until the desired 
accuracy is achieved. If the optimized sensor configuration does yield sufficiently accurate 
results, then the number of sensors used in the configuration is decreased until no more 
sensors can be removed without causing the configuration to fail the accuracy test. The 
resulting number of sensors and locations of these sensors are considered to then be the best 
estimate of the minimum number of sensors which will achieve the desired average accuracy 
averaged over all selected locations. 

Korbicz's algorithm formulation does not contain any existence theorem showing that 
the approximations so arrived at converge to a solution which is globally minimal. The 
filtering portion of the algorithm relies on a variation of the Kalman Filter theorems, but the 
sensor location portion of the algorithm is a derivation using control theory principles, and 
is demonstrated in an example in Korbicz's papers. 

5.8 Suboptimal filtering and sensor location results 

Using the previously specified contaminant concentration scenario, variations of sensor 
configurations containing sensors placed at spatial locations xi» X6, x\ \ were then applied 
and analyzed using the suboptimal filtering and sensor location algorithm. A step input was 
introduced into the simulations at spatial location xn at time 0.01 sec. The simulations 
were run for simulated durations of 2 seconds, which was enough time for the states to 
progress to near-equilibrium (the spatial locations were close together). Testing was 
performed to establish the appropriate simulation time step and other algorithm parameters 
which can be varied. The time step increment was started at the stability threshold and 
reduced until the simulation results did not appreciably change; the time step finally selected 
was 0.003 sec. The other parameter selections and variations are described in the next 
section. 
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The suboptimal algorithm is designed to assess the state prediction error of a number of 
sensors in a particular configuration under investigation, perform comparisons between 
configurations of sensors, and place those sensors in optimal locations. When three sensors 
for three locations are being assessed, only one configuration is possible. When the use of 
two sensors is being assessed for three possible sensor locations, three different sensor 
configurations are possible, and the algorithm should indicate the best two locations for 
those sensors, as will be shown below. 

Sensor configurations were tested with sensors placed at each possible location and 
then placed in groups of two. A full configuration consisted of the placement of sensors at 
each of the three possible spatial locations (*l, *6, *ll)- 1° the terminology of switching 
sensors on and off (1 or 0), this full configuration would be X = [1 1 1], since a sensor is 
placed at each possible spatial location which can accept a sensor. The placement of 
sensors only at locations x\ and x6 translates to A, = [1 10], For each sensor configuration, 
correlation neighborhoods of An = 0, 1, 2, 3, 4, 5, 11 were run. The average state error 
was calculated for each run. A summary of the results for An = 4 is shown in Table 2. 
Plots of the actual system states, sensor measurements, and the estimated system states for 
An = 4 and sensors at locations xe and x\\, are shown in Figure 4. 

Table 2. Results of Korbicz Suboptimal Filter and Sensor Location Algorithm 


Sensor 

! Correlation i 

Average 

Prioritized 

Placement 

| Neighbors | 

State 

Sensor 

A = 

I 

* 

>> 

0\ 

1 

An j 

Error 

Locations 

[1 1 11 

1 4 

0.01108 

IUJLLL 

[1 i oj 

j 4 l 

0.04481 

LteJlJl 

[i°il 

1 4 ! 

0.01460 

i hi 6 1] 

[Oil] 

I 4 1 

0.01103 

1 HI 6 11 
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Figure 4. Korbicz Suboptimal Filter & Estimate of States, Sensors at x$ and x\ j, 
Response to a Step Input at Location xil. Aw = 4 

When sensors were placed at locations xj$ and x\ i, and using Aw = 4, the average state 
error over the 2.0 sec. duration of the simulation was calculated to be 0.0110, slightly less 
than that calculated for the Extended Kalman Filter using the same sensor configuration 
(previously 0.01 199). The prioritized sensor locations calculated by the suboptimal sensor 
placement algorithm was [11 6 1], which means that the algorithm indicates that the best 
locations for two sensors for this step input scenario is in locations xi ] & X6. As shown by 
the results from the Extended Kalman Filter tests, these are probably the best locations. 
Placing sensors at spatial locations xi & xi i yielded a significantly greater average state 
error of 0.0146. 

If the simulation started off placing sensors at locations x\ & x<>, using JEn = 4, the 
resulting sensor location priority would indicate that placing sensors at X6 & X) 1 would 
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yield better results. Even though the pair of locations x\ & x\ 1 yields better estimates of 
the concentrations when using Korbicz's Suboptimal Filtering algorithm, the pair of 
locations *6 & x\ \ are the ones that actually yield the best concentration estimation results, 
since the Extended Kalman Filter is the likely Filter to be used in the space habitat. 

5.9 Sensor location algorithm discussion 

There is a large variation in the algorithm selected prioritized sensors for correlation 
neighborhoods An = 0 - 32. There are two main reasons for this, corresponding to two 
mistakes in the formulation of the adjoint equations. The first mistake is in the sign of the 
adjoint equations. The equations to calculate an adjoint variables, as derived in (W.F. 
Ramirez, 1994), are, 



(5.49) 


(5.50) 


(5.51) 

The corresponding Korbicz algorithm equations, as presented 

in Eqns. (42), (43), & 


(44), are of the opposite signs. Korbicz seems to have taken the calculations of the 
continuous-time adjoint variables from (J. Korbicz et al., 1988), and to have simply called 
them discrete-time (J. Korbicz, 1991). This sign reversal has two primary effects in the 
calculation of the adjoint variables: 1) the adjoint equations have significant terms which 
oscillate between negative and positive values, and 2) the Hamiltonian has to be maximized. 
The first point is easy to see. For example, in Eqn. (43), if we assume that y(V, k| k) is a 
very small value, and that the (N X )J term is positive-definite such as in our simulation, then 
Eqn. (43) simply changes the sign of \|f(x,k|k). This also occurs for the calculation of 
k| k) . 


57 



This sign reversal also forces Korbicz to try to maximize the Hamiltonian, when it 
should actually be minimized because we want to minimize the cost functional. At one 
iteration, the values of the X-dependent portion of the Hamiltonian could indicate that the 
sensor location priority is one particular order, and due to the sign reversals, the very next 
iteration will indicate that the sensor location priority is of the opposite order. The values 
of £j bounce between extremes as the adjoints are calculated. 

The second mistake is in the calculation of the adjoint equations. It seems that Korbicz 
was not rigorous in his calculations; the use of the (•) placeholder in the denominators of 
Eqn. (42), (43), & (44) ignores some complex subscript problems. In examining Eqns. 
(37), (39), & (33), we see that y is a function of G, which. is a function of M (different time 
index from Af(»)), which is a function of y . Since y is therefore a function of y down the 
line, then it would follow that there should be a <|>(x, k| k) in Eqn. (43) since the Hamiltonian 
H contains <|>(x, k| k) x y . 

In addition, the calculation of 0(x,x^k|k) seems to be oversimplified, if Korbicz's 
equation is taken literally (it can result in 9(x,x^k|k) being diagonal, with all of the terms 
being the same). It is therefore surprising that the calculation of the correct sensor location 
priority actually returns correct values, though only for large correlation neighborhoods. 

5.10 Conclusions 

Once the parameters are set correctly, the Korbicz Suboptimal Filter estimates the 
concentrations of the contaminants (states) as well as the Extended Kalman Filter, for the 
contamination scenario tested in this paper. 

The Sensor Location Algorithm does not work properly until the correlation 
neighborhood reaches some value, though this value is not presently able to be determined 
beforehand. Once that value is achieved, though, the Suboptimal Filter portion may imply 
that one particular sensor configuration is optimal (through examining the average state 
error), while the Sensor Location Algorithm may select a different sensor configuration as 
optimal (which coincided with the Extended Kalman Filter results in our example). 
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The calculation of the adjoint equations in the Sensor Location Algorithm is suspect; 
some terms may have been left out or approximated, though this is not mentioned by 
Korbicz. It was necessary to re-derive all of Korbicz's equations to determine the correct 
equation multiplications, since the subscripted equations presented by Korbicz in his papers 
are extremely confusing. When this re-derivation was performed, some neglected terms 
were discovered. Before the Algorithm can be used with more confidence, the significance 
of those missing terms needs to be evaluated. Until then, use this algorithm with caution. 

The general methodology developed by Korbicz looks promising, but the sensor 
placement algorithm equations presented by Korbicz are not rigorous, and they leave a 
significant amount of variation in their application. We shall be recalculating the algorithm 
equations in an attempt to obtain more consistant results. 
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Chapter 6 


Conclusions 


To date we have developed two and three dimensional distributed parameter models 
of contaminant transport, developed a new Implicit Kalman Filtering approach for 
contaminant identification, and developed a suboptimal sensor placement algorithm. 
This coming year we plan to work on a contaminant source diagnosis problem, develop 
a three dimensional contaminant visualization program, use the Implicit Kalman 
Filter to estimate three dimensional contaminant concentration profile, and develop 
an optimal sensor placement algorithm. 
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